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1  Introduction 


Fourier  techniques  have  been  a  popular  analytical  tool  in  physics  and  engineering  for  more 
than  two  centuries.  A  reason  for  this  popularity  is  that  the  trigonometric  functions  e’"*  are 
eigenfunctions  of  the  differentiation  operator  and  thus  form  a  natural  basis  for  representing 
solutions  of  many  classes  of  differential  equations. 

More  recently,  the  arrival  of  digital  computers  and  the  development  of  the  fast  Fourier 
transform  (FFT)  algorithm  in  the  1960s  (see  [7])  have  established  Fourier  analysis  as  a  powerful 
and  practical  numerical  tool.  The  FFT,  which  computes  discrete  Fourier  transforms  (DFTs), 
is  now  a  central  component  in  many  scientific  and  engineering  applications,  most  notably  in 
the  areas  of  spectral  analysis  and  signal  processing.  Numerous  applications,  however,  involve 
unevenly  spaced  data,  whereas  the  FFT  requires  input  data  to  be  tabulated  on  a  uniform  grid. 
In  this  paper,  we  present  a  collection  of  algorithms  which  overcome  this  limitation  of  the  FFT 
while  preserving  its  computational  efficiency.  These  algorithms  are  designed  for  the  efficient 
computation  of  certcun  generalizations  of  the  DFT,  namely  the  application  and  inversion  of  the 
transformation  F  :  —*  defined  by  the  formulae 

N/2-I 

Fia)j=  (1) 

k=-N/2 

for  j  —  l,...,iV,  where  x  =  {ii,.,.,z/v}  is  a  sequence  of  real  numbers  in  [— 5r,?r]  and  a  = 
{a_;v/2)  •  •  is  a  sequence  of  complex  numbers.  The  number  of  arithmetic  operations 

required  by  each  of  the  algorithms  of  this  paper  is  proportional  to 

A-logAT  +  iV-logQ^  (2) 

where  e  is  the  desired  accuracy,  compared  with  0{N^)  operations  required  for  the  direct  appli¬ 
cation  and  0{N^)  for  the  direct  inversion  of  the  transformation  described  by  (1). 

Remark  1.1  The  DFT  in  “unaliased”  form  is  described  by  the  formula 

N/2-l 

fi=  E  (3) 

k=-N/2 

for  j  =  -Nf2, . .  .,Nf2  -  1,  which  is  clearly  a  special  case  of  (1).  The  FFT  algorithm  employs 
a  sequence  of  algebraic  manipulations  to  reduce  the  number  of  operations  for  the  DFT  from 
O(iV^)  to  0{N  -log  N).  In  the  more  general  case  of  (1),  the  structure  of  the  linear  transformation 
can  also  be  exploited  via  a  combination  of  certain  analytical  results  and  the  FFT. 

The  algorithms  of  this  paper  utilize  the  fact  that  a  Fourier  series  is  a  trigonometric  polyno¬ 
mial.  When  dealing  with  the  values  of  this  polynomial  at  equispaced  nodes  on  the  unit  circle, 
the  FFT  can  be  applied.  In  this  paper,  however,  we  are  interested  in  the  values  at  nonuniformly 
spaced  nodes,  which  are  the  values  of  the  polynomial  which  interpolates  the  equispaced  values. 
The  algorithms  we  will  describe  rely  for  their  efficiency  on  a  combination  of  the  FFT  with  a  fast 
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algorithm  for  evaluating  trigonometric  polynomial  interpolants  which  uses  a  version  of  the  Fast 
Multipole  Method  (FMM)  specifically  designed  for  the  geometry  of  the  circle.  This  interpola¬ 
tion  algorithm  is  closely  related  to  the  interpolation  algorithm  described  in  [9]  for  polynomials 
tabulated  on  the  line. 

Remark  1.2  Throughout  this  paper  we  will  be  using  the  well  known  Lagrange  representation 
of  polynomial  interpolants.  For  a  function  /  :  C  — ►  C  tabulated  at  nodes  this  is 

defined  by  the  formula 

=  (-I) 

j=i  it=i 

Following  is  a  plan  of  this  paper.  Section  2  contains  a  number  of  results  from  analysis 
and  approximation  theory,  and  in  Section  3  we  describe  both  formally  and  informally  how 
these  results  are  used,  together  with  the  FMM,  in  the  construction  of  the  fast  algorithms  of  this 
paper.  Results  of  several  of  our  numerical  experiments  are  presented  in  Section  4  to  demonstrate 
the  performance  of  these  algorithms,  and  finally  in  Section  5  we  discuss  several  generalizations 
and  conclusions. 

Remark  1.3  An  alternative  approach  to  the  problems  of  this  paper  is  presented  in  [10],  where 
an  interpolation  scheme  based  on  the  Fourier  analysis  of  the  Gaussian  beU  is  used  in  place  of 
the  FMM-based  interpolation  scheme  of  the  present  paper.  We  compare  the  two  approaches  in 
Section  5. 


2  Mathematical  and  Numerical  Preliminaries 

This  section  is  divided  into  two  parts.  In  Subsection  2.1  we  present  several  identities  which  are 
employed  in  the  development  of  the  fast  algorithms  of  this  paper.  Subsection  2.2  contains  a 
collection  of  error  bounds  which  allow  us  to  perform  calculations  to  any  prescribed  accuracy. 

2.1  Analytical  Tools 

The  main  results  of  this  subsection  are  Theorems  2.3  and  2.4  which  describe  linear  transforma¬ 
tions  connecting  the  values  of  a  Fourier  series  at  two  distinct  sets  of  points.  Lemmas  2.1  and 

2.2  provide  intermediate  results  which  are  used  in  the  proofs  of  these  theorems. 

Lemma  2.1  Let  {xj, . . . ,i;v}  {yi,...,yAf}  be  sequences  of  real  numbers  on  the  interval 

[— jr,jr],  and  let  {u»i, . . . ,  iua;}  and  {zi,...,z/v}  be  sequences  of  complex  numbers  defined  by  the 
formulae 


(5) 

(6) 
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for  j  =  1,. .  .,N .  Then, 


n("’'  ~  ■**)  =  •  n  '  2*  •  sin((x/  -  yfc)/2) 


fori  =  and 


n(^j  -  ^k)  =  •  n  •  2*  •  sin((yj  -  yjk)/2) 


/or  j  =  1,...,N. 


Proof.  A  sequence  of  simple  algebraic  manipulations  and  trigonometric  identities  yields 
N  N 

Yiiwi-zk)  = 


_  gt(N-l)xi/2  .  jj  g«»*/2 . 2,- .  8in((l;  -  yk)l2) 


(N-i)/2  .  TT  1/2  ,  o, . 


JJ  zl'  ■  2i  •  sin((x/  -  yk)/2). 


Substituting  Zj  for  wt  and  yj  for  x;  in  (9),  we  also  obtain 

JJ(^j  -  ^k)  =  •  JJ  .  2i  •  sin((yj  -  yk)/2). 

fc=i  fc=i 

k^j  k^j 


The  following  lemma  describes  an  alternative  representation  of  the  well  known  Lagramge 
interpolation  formula  for  polynomials  in  the  case  when  the  interpolation  points  lie  on  the  circle. 

Lemma  2.2  Let  {xi,...,x;v}  ond  {yi,...,y7v}  ^  sequences  of  real  numbers  in  the  interval 
[— TT,  tt],  and  let  {/j, . . . ,  f^}  be  a  sequence  of  complex  numbers.  Further,  let  {lOi, . . . ,  Wff}  ond 
{zi,. .  .,zn]  be  sequences  of  complex  numbers  defined  by  the  formulae 

Wj  =  e’^  (11) 

z.  =  e’^  (12) 


for  j  =  Then, 


where  {c/}  and  {dj}  are  defined  by  the  formulae 


c;  =  JJ  sin((x/  -  yk)/2), 
k=l 

dj  =  fr  — = — 

s‘n((yi  -  yk)/2) 


for  j,l  = 

Proof.  Dividing  (7)  by  (8)  we  obtain 


N 

rr  Wl  —  Zk  _ 

2,  -  Zu 
k=l  ^  * 

Mi 

(yv_i)/2  N 

Jt=l 

g-i(x(-Vj)/2 

sta((*i  -  »,)/2) '  Mivi  -  n)m ' 

and  the  combination  of  (16)  with  the  fact  that 

^  cos((x(  -  yj)/2)  +  tsin((x/  -  yj)/2} 
sin((i,  -  yj)/2)  sin((i/  -  yj)/2) 

1 _ . 

tan((i,  -  yj)/2)  *’ 


gives  us 


S'-  ,5“  ■  5''  ■.-■"■'■■(isshss-')’  " 

where  {cj}  and  {dj}  are  defined  by  (14)  and  (15).  □ 


.  ~  ^k 

jfc=l  •'  * 


The  following  theorem  provides  a  formula  for  determining  the  values  of  a  Fourier  series  at 
a  set  of  points  in  terms  of  the  values  of  this  series  at  another  set  of  points. 

Theorem  2.3  Let  {xj, . .  .jXat}  and  {yi,-  ■  ■,yff}  be  sequences  of  real  numbers  in  the  interval 
[-7r,7r],  and  let  {a_!si/2^  •  •  •  7«n/2-i}  ^  «  sequence  of  complex  numbers.  Further,  let  {tui, . . .,  u>a^}, 
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{zi, . .  .,zn},  and  {gi,...,gff}  be  sequences  of  complex  numbers  defined  by  the  for¬ 

mulae 


Wj  =  e'^^ 

(19) 

II 

(20) 

A//2-1 

fi  =  £ 

(21) 

k=-N/2 

N/2-1 

9i  =  Oke'^^> 

(22) 

k=-N/2 

for  j  =  1 , . . . ,  jV .  Then, 


91  =  Cl- 


(tan((i,- yj)/2)  *)  ’ 


(23) 


where  {q}  are  defined  by  (14)  and  {dj}  are  defined  by  (15). 

Proof.  Let  the  polynomial  Pa  be  defined  by  the  formula 

N-l 

Pc{2)  =  ^  ak-N/2  • 

fcsO 


(24) 


The  Lagrange  interpolation  formula  relates  the  values  of  Pa  at  the  points  {in/}  to  the  values  at 
the  points  {z^}  via  the  expressions 


N 


N 


p„(u;/)= j;p,(^j).n 

j=i 


Wl  -  Zk 


Zi  —  Zk 
k=l  ^  * 


(25) 


for  /  =  1, . . iV,  and  applying  Lemma  2.2  to  (25)  we  obtain 

= u,r  ■  c,  ■  g  .  d,  ■  -  i) 

From  the  combination  of  (24)  and  (19)-(22)  we  see  that 

N/2-1 

Pa(Zj  )  =  z^'"^  -  Yj  Ok  ■  Zj  =  Z^^'^  -  fj 

k=-Nf2 


(26) 


and 


N/2~l 

Pa{wi)  =  ^  Qk-wj‘  =  ■  gi, 

k=-Nl2 


(27) 


(28) 
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and  finally  substituting  (27)  and  (28)  into  (26)  we  obtain 


(29) 


for/= 


In  the  case  when  the  points  {j/j}  are  equispaced  in  [— ir,ir],  the  interpolation  formula  of 
Theorem  2.3  has  a  simpler  form,  which  is  described  in  the  following  theorem.  The  result  of  this 
theorem  can  be  found  in  a  slightly  different  form  in  [11]. 

Theorem  2.4  Let  {arj, . .  .,1/^/}  be  a  sequence  of  real  numbers  on  the  interval  [— 7r,7r]  and  let 
{ci-N/2i  ^N/2-i}  ^  ^  sequence  of  complex  numbers.  Further,  let  {j/i , . . . ,  yiq}  be  a  sequence 
of  real  numbers  defined  by  the  formulae 


Vi  =  (i  -  1  -  Nl2)irlN 


(30) 


for  i  =  and  let  {wu---,wn},  {zi,...,2s},  and  be  se¬ 

quences  of  complex  numbers  defined  by  the  formulae 


Wi  = 


for  j  =  1,. .  .,N .  Then, 


.  fNxi\ 
gi  =  sin  {-—) 


Zj  = 

e‘S'^ 

N/2-1 

II 

k=-N/2 

N/2-1 

9}  = 

k=-N/2 

N 

r. 

(-ly  i 

ym 


-) 


(31) 

(32) 

(33) 

(34) 


(35) 


Proof.  From  the  combination  of  (30)  and  (33),  we  see  that  the  sequence  {q!_;v/2)  •  •  •■>aiNl2-\) 
is  the  discrete  Fourier  transform  of  the  sequence  {/i,. .  .,/n}.  In  other  words. 


1  ^ 


ikyj 


(36) 


i=i 


for  k  =  —N/2, . . . ,  N/2  —  1.  Let  us  now  define  the  function  /  :  [— tt,  x]  — »  C  by  the  formula 


N/2-1 

/(i)=  j: 

k=-N/2 


(37) 
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Substituting  (36)  into  (37)  and  changing  the  order  of  summation,  we  obtain 

N/2-1  .  N 

m  =  E  jv'E 

k=-N/2  i=l 
N  1  N/2-1 

=  E 

j=l  k=-N/2 

Observing  that  the  second  sum  in  the  expression  (38)  is  a  geometric  series  we  have 

g-.W(*-»,)/2  _  .•N(x-y,)/2 

^  1  _  e*(x-yj) 

k=-N/2 

smiNjx  -  yj)/2) 


e'i^-yj)/^  •  sin((x  —  yj)/2) 

=  (sin(JV(i  -  yj)/2))  •  (cot((i  -  yj)/2)  -  i) 

for  any  x  €  [-7r,7r].  The  definition  of  {yj}  now  yields 

sin(Ar(a;  -  yj)/2)  =  sm{Nx/2)co^Nyj/2)  -  cos(Nx/2)siTi{Nyj/2) 

=  sin(A^i/2)  •  (— 1)-', 

and  finally,  using  the  fact  that  gi  =  f{xi)  and  combining  (38),  (39)  and  (40)  we  obtain 


.  {Nx,\  (-1)^  ( 


tan((x/  -  yj)/2) 


-)■ 


(38) 


(39) 


(40) 

(41) 
□ 


2.2  Relevant  Facts  from  Approximation  Theory 

The  algorithms  of  this  paper  are  based  on  several  results  from  the  Chebyshev  approximation 
theory  of  the  function  l/tan(i).  These  results  are  contained  in  the  lemmas  and  theorems  of 
this  subsection,  numbered  2.8-2.14.  Analogs  of  these  results  for  the  function  l/i  can  be  found 
in  [9]. 

The  main  results  of  this  section  fall  into  two  categories.  Theorems  2.11  and  2.14  describe 
how  the  function  l/tan(x)  can  be  approximated  on  different  regions  of  the  interval 
using  Chebyshev  expansions.  Theorems  2.15,  2.16  and  2.17  provide  three  ways  of  manipulating 
these  expansions  which  are  needed  by  the  fast  algorithms  of  this  paper. 

We  begin  with  three  classical  definitions  which  can  be  found,  for  example,  in  [13],  [17]. 

Definition  2.1  The  n-th  degree  Chebyshev  polynomial  Tn{x)  is  defined  by  the  following  equiv¬ 
alent  formulae: 

Tn(x)  = 


cos(n  arccos  x) 
1 


(42) 


Definition  2.2  The  roots  of  the  n-th  degree  Chebyshev  polynomial  Tn  lie  in  the  in¬ 

terval  [—1,1]  and  are  defined  by  the  formulae 

for  k  =  I, . .  .,n.  They  are  referred  to  as  Chebyshev  nodes  of  order  n. 

Definition  2.3  We  will  define  the  polynomials  ui , . . . ,  ti„  of  order  n  —  1  by  the  formulae 

“<(')  =  ft 

fc=l  ^  * 

ki^i 

for  j  =  1, . .  .,n,  where  tk  are  defined  by  (44)- 

For  a  function  /  :  [—1, 1]  —*  C,  order  n  -  1  Chebyshev  approximation  to  /  on  the  interval 
[—1,1]  is  defined  as  the  unique  polynomial  of  order  n  —  1  which  agrees  with  /  at  the  nodes 
ti, . .  .,t„.  There  exist  several  standard  representations  for  this  polynomial,  and  the  one  we  will 
use  in  this  paper  is  given  by  the  expression 

(46) 

3=1 

For  the  purposes  of  this  paper,  Chebyshev  expansions  for  any  function  wUl  be  characterized  by 
values  of  this  function  tabulated  at  Chebyshev  nodes. 

Lemmas  2. 5-2. 7  provide  estimates  involving  Chebyshev  expansions  which  are  used  in  the 
remainder  of  this  section.  The  proof  of  Lemma  2.5  is  obvious  from  (42). 

Lemma  2.5  Let  Tn{x)  be  the  Chebyshev  polynomial  of  degree  n.  Then, 

|Tn(x)!  <  1  (47) 

for  any  x  e  [-1, 1]. 

Lemma  2.6  Let  Tn{x)  be  the  Chebyshev  polynomial  of  degree  n.  Then, 

|T„(x)l  >  I  •  |y  I"  (48) 

for  any  x  such  that  |x|  >  3. 

Proof.  From  Definition  2.1,  we  have 


|rn(x)l  = 
> 
> 

for  any  x  such  that  |x|  >  3. 


1  •  |(x  +  -  1)”  +  (x  -  y/x^  -  l)"j 

^  ■  |x  +  v'x^  -  (x/3)T  =  ^  •  |x  •  (1  +  v^)r 
1  ^  " 

2  ■  3 


(49) 

□ 
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Lemma  2.7  Let  Uj(i)  be  defined  by  (45).  Then,  for  any  x  €  [-1, 1], 


|u_j(x)|  <  1. 


(50) 


Proof.  It  is  obvious  from  (45)  that  nfitj)  =  1,  and  that  Uj{tic)  =  0  when  k  j.  In  addition, 
the  expression 

(51) 


k=i 


is  also  equal  to  1  at  tj  and  equal  to  0  at  all  other  t*.  Since  both  Uj  and  (51)  are  polynomials 
of  order  n  —  1 ,  we  have 

=  (52) 


*=1 


for  j  =  1, . . n.  Furthermore,  due  to  the  combination  of  (52)  and  the  triangle  inequality,  we 
obtain 

1^ 


|uj(x)|  = 

for  any  i  6  [-1, 1]. 

The  next  lemma  is  obvious. 


fc=i 


k=l 


(53) 

□ 


Lemma  2.8  For  any  a  €  [0,  |], 

tan  3a  >  3  ■  tan  a.  (54) 

The  following  two  lemmas  provide  preliminary  results  which  are  used  in  the  proof  of  Theo¬ 
rem  2.11. 


Lemma  2.9  Suppose  that  n  >2,  and  that  b>  0  and  xq  are  real  numbers  with  |io|  >  36.  Then, 
for  any  x, 


1  +  xxq  -  {x  -  xo)  •  ^ 
i=i 


1  -I-  btjXg 
btj  -  Xo 


Tn{x/b) 

TnixolbY 


(55) 


Proof.  Let  Q{x)  by  the  polynomial  of  degree  n  defined  by  the  formula 

Q(x)  =  1  -I-  xxo  -  (x  -  lo)  •  ^  hr  ^-xo  '  (?)  ■ 

It  follows  from  the  combination  of  (45)  and  (56)  that 

Q{btk)  =  1  +  btkXo  -  {btk  -xo)-J2  ■  “i(^*) 

=  1  -I-  64*0  -  (64  -  Xo)  •  =  0,  (57) 

04  -  Xq 
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for  ^  =  1, . . . ,  n.  Clearly,  then,  Q{x)  satisfies  the  conditions 

Qixo)  =  1  + 

Q{bti)  =  0 


It  is  clear  that  the  function 


Qibtn)  =  0. 


n  ,  ,2,  T.(x/b) 


is  also  a  polynomial  of  degree  n  which  satisfies  the  n  +  1  conditions  (58).  Therefore, 

and  (55)  follows  as  an  immediate  consequence  of  (56)  and  (60). 


Lemma  2.10  Suppose  that  n  >  2,  and  that  b  >  0  and  xq  are  real  numbers  with  jarol  ^  3^- 
Then, 


1  +  xiQ  ^  1  +  btjXp  / x\  1  -I-  96^ 

X  ~  Xq  btj  -  Xq  \b  J  ^  6-5" 


for  any  x  £  [-6, 6]. 


Proof.  Dividing  (55)  by  {x  -  lo)  and  taking  absolute  values,  we  obtain 


X  -  lo 


1  -i-  IXq  ^  ^  +  btjXQ  y  /f')  1  +  X§  |rn(x/fr)| 

I  -  iQ  btj  -  Xq  ^^\b)  ~|i-io|  |T„(io/5)r 


Due  to  Lemmas  2.5  and  2.6  we  have 


|7;(x/6)|  <  1 


for  any  i  6  [-5,6],  and 


1  +  x^o 

Tn{xo/b) 


<(l+*S)-2-|^r<i(l+(36)>) 


for  any  |*o|  >  36.  Finally,  substituting  (63)  and  (64)  into  (62),  we  obtain 


1  -f  xxq  ^  1  +  btjXp  /x\  1  -I-  96^ 

X  -  Xo  btj  -  Xq  \  6  y  ^  6-5" 


for  any  x  6  [-6,6]. 
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Theorem  2.11  Suppose  that  n  >  2,  and  that  a  and  $o  are  real  numbers  xoith  0  <  a  <  ^  and 
3a  <  l^ol  <  f  •  Then, 


^  1  +  tj  tan  a  tan  Bq  /  tan  0  \ 
2^  #  tan/i—fanflQ  ^^-^Vtana/ 


j  tar  (0  -  ^  tj  tan  a  —  tan  9^ 

for  any  9  €  [-a, a]. 


1  +  9  tan^  a 
tan  a  ■  5” 


(66) 


Proof.  Let  9  €  [—a, a].  Then,  defining  the  real  numbers  b,  x,  and  xq  by  the  formulae  6  =  tana, 
X  =  tan^  and  xq  =  tan^oj  we  observe  that  |x|  <  6,  and,  due  to  Lemma  2.8,  ixo|  >  tan  3a  >  36. 
We  also  observe  that 

1  1  +  tan  9  tan  Bq  1  +  xxq 


tan(fl  —  ^o)  tan  9  —  tan 


Xo 


and 


^  1  +  tj  tan  a  tan  Bq 


Bq  /tand\  _  ^  1  +  tj6xo  f  x\ 
'o  Vtana/  "  tjb-xo  ^^\b) 


“  tj  tana  —  tan^o 

It  follows  from  the  combination  of  equations  (67)  and  (68)  and  Lemma  2.10  that 


(67) 

(68) 


tan(fl 


1  ^  1  +  tj  tan  a  tan  /  tan  9  \ 

'  ~  ^o)  ^  tj  tan  a  — tan  00  ^  Vtana/ 


1  +  xxo  ^l  +  tj6xo  f  x\ 

- ETrrrr-iUj 


X  -  Xo  ^  tj6  -  Xo 

1  +  96^  1  +  9  tan^  a 


(69) 


6-5” 


tan  a  • 5” 


for  any  9  €  [-a,  a] 


The  following  two  lemmas  provide  preliminary  results  which  are  used  in  the  proof  of  Theo¬ 
rem  2.14. 

Lemma  2.12  Suppose  that  n  >  2,  and  that  6  >  0  and  Xq  are  real  numbers  with  |xo|  <  6.  Then, 
for  any  x. 


X  +  36x„  -  (34  -  xx„) .  ±  .  „,(x)  =  (I  +  34x„)  • 

Proof.  Let  Q(x)  be  the  polynomial  of  degree  n  defined  by  the  formula 

Q(x)  =  X  +  36xo  -  (36  -  xxq)  •  •  «;(*)• 

^  36  -  tjXo 


(70) 


(71) 
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It  follows  from  the  combination  of  (45)  and  (71)  that 

Qih)  =  +  36zo  -  (36  -  tkXo)  •  ^  ^  •  ^i{tk) 

"  3o  -  ijio 

=  tfc  +  36*0  -  (36  -  IfcXo)  •  - — -  =  0, 

do  —  IfcXo 

for  k  =  1, . .  .,n.  Clearly,  then,  Q(x)  satisfies  the  conditions 

Q(36/xo)  =  36/xo  +  36xo 
Q(ti)  =  0 


Q(tn)  =  0. 


It  is  clear  that  the  function 


Vxo  /  Tn 


Tn(x) 


»(36/xo) 

is  also  a  polynomial  of  degree  n  which  satisfies  the  n  +  1  conditions  (73).  Therefore, 

and  (70)  follows  as  an  immediate  consequence  of  (71)  and  (75). 


(72) 


(73) 

(74) 

(75) 
O 


Lemma  2.13  Suppose  that  n  >  2,  and  that  6  >  0  and  xq  ore  real  numbers  with  |xo|  <  6.  Then, 

3(1  +  62) 


X  +  36xo  ^  +  36xo  f  . 


6-5" 


1 36  -  xxq  ^  36  -  tjXQ 

for  any  x  €  [-1, 1]. 

Proof.  Dividing  (70)  by  (36-  xxq)  and  taking  absolute  values,  we  obtain 


X  +  36xo  ^  tj  +  36xo  -  ^ 

- 1-  ir—^:  ■  «^(2:) 


1 


|36-  xxol 


1 36  -  xxo  ^  36  -  tjXQ 
In  addition,  due  to  Lemmas  2.5  and  2.6  we  have 

|Tn(x)|  <  1 

for  any  x  €  [-1,1],  and 


o.. 

—  +  36xo 


\Tn{x)\ 

|r„(36/xo)|' 


36/xo  +  36xo 


Tn(36/xo) 


<  36  •  (1/xo  +  xo)  •  2 


(76) 


(77) 


(78) 


(79) 
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for  |io|  <  b.  Substituting  (78)  and  (79)  into  (77),  we  obtain 


I  +  36xo 
36  —  xxq 


^ij  +  36io 
"  36  -  tjXo 


•  Uj(x) 


3(1  +  6^) 
6-5" 


(80) 


for  any  x  [-1, 1]. 


□ 


Theorem  2.14  Suppose  that  n  >  2,  and  that  a  and  Bq  are  real  numbers  with  0  <  a  <  ^  and 
l®o|  <  a-  Then, 


1  ^  tj  +  3tanatan0o 

tan(d  -  Bo)  ^  3  tan  a  —  tj  tan 


< 


6 

sin  2a  •  5” 


(81) 


for  any  B  such  that  3a  <  |d|  <  ^. 


Proof.  Let  B  be  any  real  number  such  that  3a  <  1^1  <  Then,  defining  the  real  numbers  6,  x 
and  xq  by  the  formulae  6  =  tana,  x  =  3  tan  a/  tand  and  xq  =  tanffo>  we  observe  that  |io|  <  6, 
and,  due  to  Lemma  2.  8,  \x\  <  1.  We  also  observe  that 


and 


1  _  1  +  tan  B  tan  _  1  +  36xo/x  _  x  -f  36xo 

tan(0  -  ^o)  tand-tan^o  36/x  — xo  36-xxo’ 

“}■  3  tan  a  tan  Bq  / 3  tan  a\  tj  -!■  36x()  ,  , 

"  3  tan  a  —  tj  tan  V  tan  B  )  ~  36  -  tjXo 


It  follows  from  the  combination  of  equations  (82)  and  (82)  and  Lemma  2.12  that 


1  ^  tj  +  3tanatan0o  /3tana\ 

tan(0  -  6>o)  ^  3  tan  a  —  tj  tan  Bq  V  tan  ^  / 

X  +  36xo  tj  +  36xo  .  , 

= - >  — - -  ti  ,•  ( X ) 

36  -  xxo  ^  36  -  tjXo  ^ 

„  1  +  6^  3  sec*  a  6 

6  •  5”  tan  a  •  5”  sin  2a  •  5” 

for  any  B  such  that  3a  <  |0|  <  |. 


(82) 

(83) 


(84) 

(85) 
□ 


The  following  three  theorems  provide  formulae  for  translating  along  the  interval  [— f»f] 
Chebyshev  expansions  of  the  type  described  in  the  previous  two  theorems.  Theorem  2.15  pro¬ 
vides  a  formula  for  translating  expansions  described  in  Theorems  2.11,  Theorem  2.16  describes 
a  mechanism  of  converting  the  expansion  of  Theorem  2.14  to  the  expansion  of  Theorem  2.11, 
and  Theorem  2.17  provides  a  way  of  translating  the  expansion  of  Theorem  2.14. 
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Theorem  2.15  Suppose  that  n,N  >2,  and  let  o,c,d  be  real  numbers  such  that  0  <  a  <  ir/8 
and  [c  —  d, c  +  d]  C  [—a,  a].  Let  the  function  /  :  [— f » f  ]  C  ^  defined  by  the  formula 


/W  =  i;s;;(fbo 


where  3a  <  |5fc|  <  ^  for  k  =  and  ai,...,aN  is  a  set  of  complex  numbers.  Further, 

let  be  a  set  of  complex  numbers  defined  by  the  formula 

9k  =  /(arctan(f/ttana))  (87) 

for  fc  =  1, . . . ,  n,  and  let  9i, . .  .,9n  be  a  set  of  complex  numbers  defined  by  the  formula 


i=i 


tan(c  +  arctan(t;t  tan  d))' 


for  fc  =  1, . . . ,  n.  Then,  for  any  6  £  [c  —  d,c  d\. 


(n  +  1)(1  +  9tan^  a) 
tan  a  •  5"  ’ 


where  A  = 

Proof.  It  follows  from  the  triangle  inequality  that 


where 


5i  =  f{9)  -  ^2  fi^  +  arctan(<jfc  tan  d))  •  Uj  ^ 


tan(^  -  c)' 


CtilU 

52  =  arctan(4  tan  d))  -  9k)  ■  «j  (~^”an  (f~)  |  ‘ 

Combining  Theorem  2.11  with  Lemma  2.7  and  the  triangle  inequality,  we  have 


5i  <>1 


1+9  tan*  a 
tan  0-5”  ’ 


n  n  y 

S2  <  ^2  /(c  +  arctan(<fc  tan  d))  -  E».«i 

fc=i  j=i  ' 


tan(c  +  arctan(/*  tand))' 


<  An 


1  +  9  tan*  a 
tan  a  •  5" 
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where  A  =  I®* I-  Finally,  substituting  (93)  and  (94)  into  (90)  we  obtain 

1  +  9  tan^  a 


for  any  0  6  [c  -  d,  c  +  d]. 


<  y4  •  (n  +  1)  ■ 


tan  a  •  5" 


(95) 


Theorem  2.16  Suppose  that  n,N  >2,  and  let  a,c^d  be  real  numbers  such  that  Q  <  a  <  t/8 
and  |c|  —  d  >  3a.  Let  the  function  f  :  [-f  ,f]  — ►  C  6e  defined  by  the  formula 


^  tan(«  -  0k) 


(96) 


where  0k  6  [—a,  a]  for  k  =  1,.  ..,N,  and  ai,...,ajv  is  a  set  of  complex  numbers.  Further,  let 
be  a  set  of  complex  numbers  defined  by  the  formula 

=  /(arctan(3tan(a)/tik))  (97) 

for  k  =  1, . . . ,  n,  and  let  9i, . .  .,9n  be  a  set  of  complex  numbers  defined  by  the  formula 

3  tan  a  \ 


9k 


j=l 


^tan(c  +  arctan(tjt  tan  d))J 
for  k  =  1, . . n.  Then,  for  any  ^  €  [c  —  d,c  +  d], 

3n  sec^  a  +  1  +  9  tan^  a 


(98) 


tan  a  •  5” 


(99) 


where  A  =  J2k=\ 

Proof.  It  follows  from  the  triangle  inequality  that 


where 


and 


5i  = 


52  = 


/(^)  -  II  /(c  +  arctan(tfc  tan  d))  •  uj  |  ’ 

53(/(c  + arctan(tfctand))  -  9k)-Uj 


Combining  Theorem  2.11  with  the  triangle  inequality  gives  us 

1  +  9  tan^  a 


Si  <  A 


tan  a  ■  5"  ’ 


(100) 

(101) 

(102) 

(103) 
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and  from  the  combination  of  Theorem  2.14,  Lemma  2.7  and  the  triangle  inequality,  we  have 


S2  <  ^ 

k=l 

<  An  • 


«  / 

/(c  +  arctan(tfc  tan  d))  - 

i=i 

3  sec^  a 


3  tan  a 


tan(c  +  arctan(tit  tan  </))> 


tan  a  •  5" ’ 

where  A  =  J2k=i  Finally,  substituting  (103)  and  (104)  into  (100)  we  obtain 

3n  sec^  a  +  1  +  9  tan^  a 


for  any  0  £  [c  -  d,c  +  d]. 


<  A 


tan  a  •  5” 


(104) 


(105) 


Theorem  2.17  Suppose  that  n,N  >2,  and  let  a,c,d  be  real  numbers  such  that  0  <  a  <  7r/8 
and  [c  -  d,c  +  d]  D  [-a, a].  Let  the  function  f  :  [— ^,  C  be  defined  by  the  formula 


N 


m  =  i: 


Qk 


“  t!in{0- 0k) 


(106) 


where  0k  6  [—a,  a]  for  k  =  1 , . . . ,  ,  and  oi , . . . ,  ojv  **  o  set  of  complex  numbers.  Further,  let 
$1 , . . . ,  $n  be  a  set  of  complex  numbers  defined  by  the  formula 

4>k  =  /(arctan(3tan(a)/tjt))  (107) 

for  k  =  1, . . . ,  n,  and  /et  $1, . . . ,  be  a  set  of  complex  numbers  defined  by  the  formula 


^k 


i=l  ^ 


3  tan  a 


A 


^tan(c  +  arctan(3  tasi{d)/tk))  J 
for  k  =  1, . . . ,  n.  Then,  for  any  0  such  that  |^  —  c]  >  3d, 

3(n  +  1)  sec^  a 


(108) 


tan  a  ■  5" 


(109) 


where  A  =  Ylk=\  lo'fel- 

Proof.  It  follows  from  the  triangle  inequality  that 


<  5i  +  52 


where 


5,  = 


/(^)  -  ]^/(c  +  arctan(3tan(d)/tit))-«j  (ta.n(g"-^c))  |  ’ 
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ana 

*^2  =  ^  (/(c  +  arctan(3  tan(d)/4))  -  ^k)  ■  «j  |  •  (112) 

Combining  Theorem  2.14  with  Lemma  2.7  and  the  triangle  inequality,  we  have 

.  3  sec*  a 


Si  <A 


tan  a  •  5”  ’ 


<  y4n  • 


3  sec*  a 
tan  a  •  5”  ’ 


where  A  =  |ofc|.  Finally,  substituting  (113)  and  (114)  into  (110)  we  obtain 


3(n  +  1)  sec*  a 
tan  a  •  5”  ’ 


(114) 


(115) 


for  any  0  such  that  \0  —  c\>  3d. 


3  Application  of  the  FMM  to  Nonequispaced  FFTs 

For  the  remainder  of  this  paper  we  will  be  considering  the  mapping  F  :  defined  by 

the  formulae 

Ar/2-i 

F{a)j=  ^  (116) 

k--N/2 

for  j  =  1,...,A^,  where  x  =  is  a  sequence  of  real  numbers  in  [— 7r,7r]  and  q  = 

{ai, . .  .,a^}  is  a  sequence  of  complex  numbers.  We  are  interested  in  the  efficient  application 
and  inversion  of  the  transformation  F  and  its  transpose.  More  formaUy,  we  will  consider  the 
following  four  problems: 

•  Problem  1:  Given  a,  find  /  =  F(a). 

•  Problem  2;  Given  a,  find  /  =  F^{a). 

•  Problem  3:  Given  /,  find  q  =  F~^{f). 

•  Problem  4:  Given  /,  find  a  =  (f^)~*(/). 

This  section  consists  of  four  parts.  In  Subsection  3.1  we  describe  briefly  how  the  one  di¬ 
mensional  Fast  Multipole  algorithm  of  [9]  can  be  applied  to  the  problems  of  this  paper,  in 
Subsection  3.2  we  outline  a  set  of  four  algorithms  for  these  problems.  Subsection  3.3  contains 
more  formal  descriptions  of  these  algorithms,  and  finally  in  Subsection  3.4  we  discuss  a  gener¬ 
alization  of  Problems  1-4. 
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3.1  FMM  and  Trigonometric  Interpolation 

There  exist  a  number  of  different  formulations  of  the  trigonometric  interpolation  problem  (see 
[17]).  The  version  we  will  use  for  the  purposes  of  this  paper  is  described  as  follows:  given  a 
set  of  points  {yj , . . . ,  and  function  values  {/i , . . . ,  /at},  evaluate  the  interpolating  Fourier 
senes  at  the  points  {xi, . .  .,ia/}.  According  to  Theorem  2.3,  these  values  are  given  by  the 
formulae 


9i  Cl  g/,  dj  *) 

for  /  =  1, . . . ,  A,  where  {c/}  and  {dj}  are  defined  by  the  formulae 

N 

c/  =  JJ  sin((x,  -  yk)f2)  =  eEir=i  M*m((ii-yfc)/2)) 


(117) 


(118) 


for  /  =  1,...,A,  and 


N  , 

d;  =  TT - = - £-Efc=i,fc^jM»>n((yj-y»)/2)) 

iiSin((y,  -yfc)/2)  " 


for  j  = 


(119) 


Remark  3.1  The  FMM  algorithms  of  [9]  are  designed  to  evaluate  expressions  of  the  form 

N 

EOfc 

in  O  ^Alog^-jj  arithmetic  operations,  where  e  is  the  desired  accuracy.  With  a  few  minor 
modifications  they  can  also  be  used  to  evaluate  expressions  of  the  form 


N 

tan(x  -  ifc) 


(121) 


^ln(sin(x  -  xjt)),  (122) 

a,nd  hence  expressions  of  the  form  (117)  for  the  same  computational  cost.  Moreover,  the  algo¬ 
rithmic  procedure  for  the  kernel  1/ tanx  is  virtually  identical  to  that  for  1/x,  and  the  various 
expansions  required  by  the  algorithms  are  manipulated  via  Theorems  2.15,  2.16  and  2.17  (see 
[9]  for  detailed  descriptions  of  these  algorithms). 
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3.2  Informal  Descriptions  of  the  Algorithms 

In  this  subsection  we  outline  how  a  fast  trigonometric  interpolation  scheme  can  be  used  to 
construct  efficient  algorithms  for  Problems  1-4  of  this  paper. 

We  begin  with  some  notation. 

T  :  will  denote  the  matrix  which  maps  a  sequence  of  N  complex  numbers  to  its 

discrete  Fourier  transform.  T  is  defined  by  the  formulae 


_  g2ir.  (i-Ar/2-l)-(fc-JV/2-l)/Af 


(123) 


for  j,  fc  =  1, . . . ,  and  it  is  well  known  that  =  T ^  and  that  =  ^!F. 

Remark  3.2  T  and  can  each  be  applied  in  O(JVlogiV)  operations  via  the  FFT. 

V  :  will  denote  the  matrix  which  maps  the  values  of  an  JV-term  Fourier  series  at 

N  equispaced  points  {yi, . . . ,  vm}  on  [— tt,  tt]  to  the  values  of  this  series  at  the  arbitrarily  spaced 
points  {xi, . .  .,x/v}.  According  to  Theorem  2.4,  V  is  defined  by  the  formulae 


^  .  /Ax.\  (-1)*  /  1  \ 

2  j'  N  \tan((x,  -  yO/2)  ‘j 

ioT  j,k  =  1, . . . ,  AT.  It  follows  directly  from  (124)  that 


(124) 


Pj  =  .  sin 


( 


(Nxk\ _ 

V  2  /  I  tan((xfc  -  yj)/2) 


-) 


(125) 


foT  j,k  =  I,..  .,N.  The  inverse  of  the  mapping  V  converts  the  values  of  an  A-term  Fourier  series 
at  the  points  {ii,...,i;v}  to  the  values  of  this  series  at  the  equispaced  points  {yi,-..,yN}- 
'P~^  is  therefore  given  analytically,  and  according  to  Theorem  2.4  it  is  defined  by  the  formulae 


'Pjk  (tan((yj  -  Xfc)/2)  ‘) 


(126) 


ior  j,k  =  1,...,A,  where  ci,...,c/v  and  d\,...,dN  are  sequences  of  real  numbers  defined  by 
the  formulae  (118)  and  (119).  It  follows  directly  from  (126)  that 


‘^^•"*'-(tan((yfc-x,)/2)  ’) 


(127) 


for  =  1,...,  A. 

Remark  3.3  P,  ^  V~^  and  are  all  of  the  same  form,  and  each  can  be  applied  with 

a  relative  precision  £  in  O  ^Alog  (7))  operations  via  the  FMM  (see  Section  3.1). 
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Observation  3.4  From  the  combination  of  (116),  Theorems  2.3  and  2.4,  and  several  elemen¬ 
tary  matrix  identities,  we  see  that 

F  =  VF 

fT  ^ 

F-^  =  (128) 

z= 

Furthermore,  due  to  Remarks  3.2  and  3.3,  F,  F^,  F~^  and  (F^)~^  can  each  be  applied  in 
O  (iVlog  arithmetic  operations. 

3.3  Formal  Descriptions  of  the  Algorithms 

Following  are  detailed  descriptions  of  the  four  algorithms  of  this  paper. 

Algorithm  1 

Step  Complexity  Description 

1  0{N  \ogN)  Comment  [Evaluate  Fourier  series  at  equispaced  points  using  FFT.] 

Compute  Qj  =  Eri-N/2  for  i  =  1, . , , ,  AT. 

2  0(A^log(i))  Comment  [Interpolate  in  space  domain  ] 

do  j  =  1 ,  N 

1/;  =  S;  •  (-l)VA^ 

end  do 

Compute  ft  =  9j/td.n((xi  —  yj)/2)  for  /  =  1, . . . ,  A*  using  FMM. 
dol=l,N 

ft  =  fl-i-  Ef=i  9j 
end  do 
do  /  =  l.AT 

ft  =  ft  ■  sin{Nxt/2) 
end  do 

Total  0(yVlogjV  +  iVlog(l)) 

Algorithm  2 

Step  Complexity  Description 

1  ^^(^log(l))  Comment  [Interpolate  in  frequency  domain.] 

do  j  =  1 ,  AT 

Qj  =  Qj  •sin(A^i;j/2) 

end  do 

Compute  at  =  -  ®>/  for  I  =  1, . . . ,  A^  using  FMM. 

do  /  =  1,N 

En 
>=1 

end  do 
do  /  =  1,  A 

at  =  at  ■  (-1)'/^ 
end  do 
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2  0{N  log  N)  Comment  [Evaluate  Fourier  series  at  equispaced  points  using  FFT.] 

Compute  fj  =  Ylk=SN/2  “t  •  for  j  =  1, ....  AT. 

Total  0{N  logN  +  N-  log(i)) 

Algorithm  3 

Step  Complexity  Description 

1  0(AAlog(j))  Comment  [Interpolate  in  space  domain.] 
do  j  —  l,N 

f:=frdi 
end  do 

Compute  ai  =  —  *j)/2)  for  /  =  1, . . . ,  TV  using  FMM. 

do/=  1,TV 

En  * 

j=l  fj 

end  do 
do  /  =  1,AT 

ai  =  ai  ■  Cl 

end  do 

2  0{N  log  N)  Comment  [Obtain  Fourier  coeffients  using  FFT.j 

Compute  Oj  =  i  a*  •  for  j  =  —TV/2, . . . ,  TV/2  —  1. 

Total  0(TV  •  log  TV  +  TV  •  log(i)) 

Algorithm  4 

Step  Complexity  Description 

1  0( TV  log  TV)  Comment  [Obtain  Fourier  coeffients  using  FFT.] 

Compute  dj  =  ji  ■  ICjLi  /*  ‘  for  j  =  1, ....  TV. 

2  0(TVlog(^))  Comment  [Interpolate  in  frequency  domain.] 

do  J  =  1 ,  TV 

end  do 

Compute  Off  =  —  aj/ t2m((zf  —  j/;)/2)  for  /  =  1, . . . ,  TV  using  FMM. 
do/=  1,TV 

EN 

j=i°j 

end  do 
do  /=  1,TV 

Of  =  Of  df 

end  do 

Total  0(TV  •  log  TV  +  TV  •  log(i)) 

3.4  FFTs  for  Complex  Data  Points 

Various  generalizations  of  the  problems  addressed  in  this  paper  are  mentioned  briefly  in  Sec¬ 
tion  5.  One  of  the  generalizations  of  Problems  1-4  merits  special  attention,  and  is  discussed  in 
this  section:  this  is  the  case  when  the  points  {ij}  are  complex,  and  lie  slightly  off  the  real  axis. 
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We  are  interested  here  in  the  transformations  described  by  the  formulae 

N/2-1 

/j  =  E  (129) 

k=-N/2 

for  j  =  1, . . . ,  iV,  which  is  a  generalization  of  (116)  with 

Xj  =  Tj  +  isj.  (130) 

Algorithms  1-4  can  be  modified  to  evaluate  expressions  of  the  form  (129),  provided  that  the  Sj 
are  small  (on  the  order  of  ^). 

Problems  of  this  type  are  frequently  encountered  in  signal  analysis,  computational  complex 
analysis  and  several  other  areas. 


4  Implementation  and  Numerical  Results 

We  have  written  FORTRAN  implementations  of  the  algorithms  of  this  paper  using  double 
precision  arithmetic,  and  have  applied  these  programs  to  a  variety  of  situations.  This  section 
contains  results  of  four  of  our  numerical  experiments  demonstrating  the  performance  of  our 
implementations  of  Algorithms  1-4. 

Several  technical  details  of  our  implementations  appear  to  be  worth  mentioning  here: 

1.  Each  implementation  consists  of  two  main  subroutines:  the  first  is  an  initialization  stage 
in  which  the  elements  of  the  various  matrices  employed  by  the  algorithms  are  precom¬ 
puted  and  stored,  and  the  second  is  an  evaluation  stage  in  which  these  matrices  are 
applied.  Successive  application  of  the  linear  transformations  to  multiple  vectors  requires 
the  initialization  to  be  performed  only  once. 

2.  The  algorithms  of  this  paper  all  require  the  evaluation  of  sums  of  the  form 

N/2-1 

E  (131) 

k=-N/2 

for  j  =  -Nf2, . .  .,N/2  —  1,  whereas  most  FFT  software  computes  sums  of  the  form 


/i  =  E  (132) 

fc=0 

for  J  =  0, . .  .,iV  -  1.  We  used  a  standard  FFT  to  evaluate  sums  of  the  form  (131)  by 
defining  a*  =  a*  for  /:  =  0, . . . ,  N/2  -  1,  a*  =  Ok-N  for  ^  =  ^/2, . .  .,N  —  I,  fj  =  fj  for 
j  =  0, . . . , N/2  -  1  and  fj  =  fk-N  for  j  =  N/2,. .  .,N  -  This  substitution  converts  the 
form  (131)  to  the  form  (132). 
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3.  The  algorithms  of  this  paper  require  an  FFT  of  size  proportional  to  N,  and  will  thus 

perform  efficiently  whenever  the  FFT  does.  This  restriction  on  problem  size  can  be 
removed  by  extending  the  input  vector  to  length  (i.e.  the  smallest  power  of  2 

which  is  greater  than  N)  and  padding  it  with  zeroes.  This  ensures  that  the  algorithms 
will  perform  efficiently  for  any  choice  of  IV.  In  our  implementations  these  changes  were 
made. 

4.  Each  of  the  algorithms  described  in  Section  3  utilizes  a  version  of  the  one  dimensional 
FMM.  The  version  used  in  our  implementations  was  Algorithm  3.2  of  [9]  which  was  chosen 
to  mzLximize  both  efficiency  and  accuracy. 

Our  implementations  of  the  algorithms  of  this  paper  have  been  tested  on  the  the  Sun 
SPARCstation  1  for  a  variety  of  input  data.  Four  experiments  are  described  in  this  section 
and  their  results  are  summarized  in  Tables  1-4.  These  tables  contain  error  estimates  and  CPU 
time  requirements  for  the  algorithms,  with  all  computations  performed  in  double  precision 
arithmetic. 

The  table  entries  are  described  below. 

•  The  first  column  in  each  table  contains  the  problem  size  N,  which  was  chosen  to  be  a 
power  of  2  ranging  from  128  to  2048  for  each  example. 

•  The  second  column  in  each  table  contains  the  relative  oo-norm  error  defined  by  the  formula 

Eoo=  max  max  \fj\  ,  (133) 

where  the  vector  /  is  the  algorithm  output  and  the  vector  /  is  the  result  of  a  direct 
cadculation. 

•  The  third  column  in  each  table  contains  the  relative  2-norm  error  defined  by  the  formula 

^2  =  (134) 

N>=i  /  j=i 

where  the  vector  /  is  the  algorithm  output  and  the  vector  /  is  the  result  of  a  direct 
calculation. 

•  The  fourth  and  fifth  columns  in  each  table  contain  CPU  timings  for  the  initialization  and 
evaluation  stages  of  the  algorithm. 

•  The  sixth  column  in  each  table  contains  CPU  timings  for  the  corresponding  direct  calcu¬ 
lation. 

•  The  last  column  in  Tables  1  and  2  contains  CPU  timings  for  an  FFT  of  the  same  size. 

Remark  4.1  Our  implementations  of  the  direct  methods  for  Examples  1  and  2  were  optimized 
by  using  the  fact  that  to  reduce  the  number  of  complex  exponential  computations. 
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Remark  4.2  Standard  LINPACK  Gaussian  Elimination  subroutines  were  used  as  the  direct 
methods  for  comparing  timings  in  Examples  3  and  4.  Estimated  timings  are  presented  for  larger 
N,  where  this  computation  became  impractical. 

Following  are  the  descriptions  of  the  experiments,  and  the  tables  of  numerical  results. 

Example  1. 

Here  we  consider  the  transformation  F  :  of  Problem  1,  defined  by  the  formula 

N/2-1 

F{a),=  52  (135) 

k=-N/2 

for  j  =  l,...,iV.  In  this  example,  {ri,...,i;v}  were  randomly  distributed  on  the  interval 
[— x,7r]  and  {oi-n/2^  •  •  •)®N/2-i}  were  complex  numbers  randomly  chosen  from  the  unit  square 

0  <  Re(2)  <  1,  0  <  Im(z)  <  1.  (136) 

The  results  of  applying  Algorithm  1  to  this  problem  are  presented  in  Table  1. 

Example  2. 

Here  we  consider  the  transformation  F^  :  of  Problem  2,  defined  by  the  formula 

N 

k=l 

for  j  =  -A/2,...,  A/2  -  1.  In  this  example,  {ii,...,iAf}  were  randomly  distributed  on  the 
interval  [-^r,  x]  and  (oi , . . . ,  o/v}  were  complex  numbers  randomly  chosen  from  the  unit  square 

0  <  Re(2)  <  1,  0  <  Im(z)  <  1.  (138) 

The  results  of  applying  Algorithm  2  to  this  problem  are  presented  in  Table  2. 


Example  3. 

Here  we  consider  the  transformation  of  Problem  3  where  F  is  defined  by  the 

formula 

tV/2-l 

F(a),  =  52  (139) 

k=-N/2 

for  j  =  1, . . . ,  A.  In  this  example,  (ij , . . . , ia/}  were  defined  by  the  forr  ulae 


Xj  =  -IT  +  2ir 


j  +  0.5  + 
A 


(140) 


for  j  =  1,...,A,  where  Sj  were  randomly  distributed  on  the  interval  [-0.1,0.!].  In  addition, 
{oi-N/2f  ■  •  ■jOiN/2-i}  were  complex  numbers  randomly  chosen  from  the  unit  square 


0  <  Re(z)  <  1,  0  <  Im(z)  <  1, 


(141) 
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and  the  numbers  {/j, . . . ,  f^}  were  computed  directly  in  double  precision  arithmetic  according 
to  the  formula  fj  =  F{a)j.  The  vector  /  was  then  used  as  input  for  Algorithm  3.  Results  of 
this  experiment  are  presented  m  Table  3. 


Example  4. 

Here  we  consider  the  transformation  -+  of  Problem  4  where  is  defined  by 

the  formula 

N 

=  (142) 

k=l 

for  j  =  —A/2, . . . ,  A/2  —  1.  In  this  example,  {xi, . . . ,xn}  were  defined  by  the  formulae 


-TT  +  2t  • 


j  +  0.5  + 
A 


(143) 


for  j  =  1,. . A,  where  6j  were  randomly  distributed  on  the  interval  [-0.1, 0.1].  In  addition, 
{oi,. .  .,qn]  were  complex  numbers  randomly  chosen  from  the  unit  square 


0  <  Re(2)  <  1,  0  <  Im(^)  <  1,  (144) 

and  the  numbers  {f-N/2^  •  ■  fN/2-\}  were  computed  directly  in  double  precision  arithmetic 

according  to  the  formula  fj  =  F^{a)j.  The  vector  /  was  then  used  as  input  for  Algorithm  4. 
Results  of  this  experiment  are  presented  in  Table  4. 


The  following  observations  can  be  made  from  Tables  1-4,  and  are  in  agreement  with  results 
of  our  more  extensive  experiments  for  this  particular  computer  architecture,  implementation 
and  range  of  A. 

1.  All  of  the  algorithms  permit  high  accuracy  to  be  attained,  and  the  observed  errors  are  in 
accordance  with  the  theoretically  obtained  error  bounds. 

2.  As  expected,  the  CPU  timings  for  aU  the  algorithms  grow  slightly  faster  than  linearly 
with  the  problem  size  A. 

3.  The  timings  for  Algorithms  1-4  are  similar,  which  is  to  be  expected  since  these  four 
algorithms  are  so  closely  related. 

4.  The  initialization  times  for  Algorithms  1  and  2  are  less  than  those  for  Algorithms  3  and 
4.  This  is  because  the  former  pair  does  not  incur  the  additional  cost  of  computing  the 
numbers  {c/t}  and  {dfc}. 

5.  The  evaluation  stages  of  Algorithms  1-4  are  about  15  times  as  ostly  as  an  FFT  of  the 
same  size. 

6.  Algorithms  1  and  2  can  compete  with  the  direct  method  at  about  A  =  32  ignoring 
initialization  time,  and  at  A  =  1024  including  the  initialization.  Algorithms  3  and  4  are 
always  dramatically  faster  than  the  direct  calculation  (15000  times  faster  at  A  =  2048) 
if  we  ignore  initialization  time,  and  break  even  with  it  at  A  =  64  if  we  include  the 
initialization. 
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N 

Errors 

Timings  (sec.) 

Eoo 

El 

Alg.  Init. 

Alg.  Eval. 

Direct 

EFT 

128 

0.379  E-14 

0.704  E-14 

1.04 

HIM 

0.09 

0.002 

256 

0.398  E-14 

0.116  E-13 

2.03 

■SB 

0.33 

0.005 

512 

0.499  E-14 

0.195  E-13 

3.26 

1.24 

0.012 

1024 

0.318  E-13 

0.625  E-13 

4.97 

■SB 

4.93 

0.026 

2048 

0.763  E-13 

0.204  E-12 

8.07 

19.62 

0.059 

Table  1;  Example  1,  Numerical  Results  for  Algorithm  1. 


N 

Errors 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

FFT 

128 

0.206  E-14 

0.800  E-14 

1.03 

0.08 

jljlffi 

256 

0.323  E-14 

0.136  E-13 

2.05 

■SB 

0.31 

512 

0.153  E-13 

0.343  E-13 

3.21 

1.20 

1024 

0.180  E-13 

0.654  E-13 

5.11 

■■I 

4.76 

2048 

0.470  E-13 

0.221  E-12 

8.16 

18.93 

Table  2;  Example  2,  Numerical  Results  for  Algorithm  2. 


7.  The  initialization  stage  is  much  more  costly  than  the  evaluation  stage  for  all  of  the  algo¬ 
rithms.  Implementing  the  algorithms  in  two  stages  thus  gives  considerable  time  savings 
whenever  the  same  linear  transformation  is  to  be  applied  to  multiple  vectors. 
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N 

Errors 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

FFT 

128 

0.117  E-13 

0.800  E-14 

1.28 

0.034 

2.96 

0.002 

256 

0.196  E-13 

0.137  E-13 

2.51 

0.082 

23.6 

0.005 

512 

0.344  E-13 

0.230  E-13 

4.33 

0.175 

189 

0.012 

1024 

0.107  E-12 

0.757  E-13 

7.45 

0.409 

1512  (est.) 

0.026 

2048 

0.357  E-12 

0.247  E-12 

12.97 

0.819 

12096  (est.) 

0.059 

Table  3:  Example  3,  Numerical  Results  for  Algorithm  3. 


N 

Errors 

Timings  (sec.) 

Eoo 

E2 

Alg.  Init. 

Alg.  Eval. 

Direct 

FFT 

128 

0.134  E-13 

0.806  E-14 

1.26 

0.033 

2.96 

0.002 

256 

0.511  E-13 

0.179  E-13 

2.47 

0.080 

23.6 

0.005 

512 

0.870  E-13 

0.373  E-13 

4.24 

0.173 

189 

0.012 

1024 

0.178  E-12 

0.811  E-13 

7.29 

0.407 

1512  (est.) 

0.026 

2048 

0.942  E-12 

0.369  E-12 

12.80 

0.820 

12096  (est.) 

0.059 

Table  4;  Example  4,  Numerical  Results  for  Algorithm  4. 

5  Conclusions  and  Generalizations 

In  this  paper  we  have  described  a  set  of  four  algorithms  for  computing  FFTs  for  nonequispaced 
data  to  any  specified  precision.  An  alternative  group  of  algorithms  for  the  problems  considered 
in  this  paper  is  presented  in  [10].  Similarities  and  differences  between  the  two  approaches  are 
summarized  below. 

1.  Both  sets  of  algorithms  use  a  standard  FFT. 

2.  Both  sets  of  algorithms  use  interpolation  formulae  to  transform  function  values  from 
equispaced  to  nonequispaced  points  and  vice-versa.  The  algorithms  of  this  paper  use  an 
interpolation  scheme  based  on  the  FMM,  while  the  algorithms  of  [10]  use  an  interpolation 
scheme  based  on  the  Fourier  analysis  of  Gaussian  bells. 

3.  For  the  application  of  the  linear  transformations  being  considered  the  algorithms  of  [10] 
are  the  more  efficient  of  the  two. 

4.  For  the  inversion  of  these  linear  transformations,  the  direct  schemes  of  this  paper  are 
generally  more  efficient  than  the  iterative  schemes  of  [10]  whose  complexity  is  dependent 
on  the  distribution  of  the  nodes. 

5.  The  FMM-based  approach  leads  to  a  set  of  closely  related  forward  amd  inverse  algorithms 
which  can  be  generalized  to  complex  data  points. 
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In  conclusion,  a  group  of  algorithms  has  been  presented  for  the  rapid  application  and  in¬ 
version  of  matrices  of  the  Fourier  kernel.  These  problems  can  be  viewed  as  generalizations  of 
the  discrete  Fourier  transform,  and  the  algorithms,  while  making  use  of  certain  simple  results 
from  analysis,  are  very  versatile,  and  have  a  broad  range  of  applications  in  many  branches  of 
mathematics,  science  and  engineering. 

The  results  of  this  paper  are  currently  being  applied  to  problems  in  a  diversity  of  areeis.  Ex¬ 
amples  include  problems  in  the  numerical  solution  of  parabolic  partial  differential  equations,  the 
analysis  of  seismic  data,  the  modelling  of  semiconductors,  weather  prediction  and  the  numerical 
simulation  of  fluid  behavior. 

Several  obvious  generalizations  of  the  results  of  this  paper  are  discussed  below. 

1.  Problems  1  and  2  involve  the  evaluation  of  an  IV-term  series  at  N  points.  Straightforward 
modifications  to  Algorithms  1  and  2  will  allow  the  efficient  evaluation  of  these  A^-term 
series  at  M  points,  where  M  ^  N .  These  modifications  have  been  implemented. 

2.  The  algorithms  of  this  paper  are  based  on  a  special  case  of  a  more  general  idea,  namely 
the  adaptive  use  of  interpolation  techniques  to  speed  up  large  scale  computations.  Other 
examples  of  this  approach  include  the  use  of  wavelets  for  the  construction  of  fast  numerical 
algorithms  (see,  for  example,  [1],  [4]),  and  the  use  of  multipole  or  Chebyshev  expansions 
for  the  compression  of  certain  classes  of  linear  operators  (see,  for  example,  [2],  [6],  [16]). 

3.  One  of  the  more  far-reaching  extensions  of  the  results  of  this  paper  is  a  set  of  algorithms 
for  discrete  Fourier  transforms  in  two  and  three  dimensions.  Detailed  investigations  into 
higher  dimensional  problems  of  this  type  are  currently  in  progress  and  will  be  reported 
at  a  later  date. 
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